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Abstract 

This  study  presents  the  development  and  evaluation  of  a 
numerical  solution  technique  used  to  analyze  the  two 
dimensional  transonic  airfoil  problem.  The  full  potential 
equation  of  motion  and  the  irrotationality  condition  are  used 
for  the  governing  equations.  A  coordinate  transformation  is 
applied  to  the  governing  equations  to  limit  the  domain  of  the 
problem  and  for  easier  application  of  the  boundary 
conditions.  The  method  of  lines  is  then  used  to  reduce  the 
equations  from  partial  to  ordinary  differential  equations. 
Solutions  obtained  using  this  numerical  solution  technique 
are  compared  to  published  data  for  both  incompressible  and 
compressible  flowfield  cases.  It  is  concluded  that  the 
solution  technique  developed  in  this  study  is  accurate  and 
efficient  when  analyzing  subcritical  flowfields  around 
circular  airfoil  shapes.  _ 


APPLICATION  OF  A  FINITE  DIFFERENCE  METHOD 


TO  THE  TRANSONIC  AIRFOIL  PROBLEM 
I.  Introduction 

Background 

The  development  of  an  accurate  and  efficient  numerical 
solution  technique  for  the  2-D  transonic  airfoil  problem  is 
one  of  the  greater  challenges  of  modern  aerodynamic  research. 
The  problem  is  physically  and  mathematically  complex  hence 
the  difficulty  in  developing  a  model  which  accurately 
reflects  what  is  physically  occurring  in  the  transonic  flow- 
field.  Nevertheless,  recent  trends  in  aircraft  design  and 
performance  requirements  have  created  the  need  for  accurate 
solutions  to  the  transonic  airfoil  problem. 

To  appreciate  the  difficulty  in  developing  a  solution 
technique  for  this  problem  one  must  have  an  understanding  of 
the  nature  of  the  physical  phenomena  involved.  In  the 
transonic  airfoil  problem  there  is  a  mixture  of  subsonic 
(M<1)  and  supersonic  ( M > 1 )  flow  regions  giving  rise  to  what 
is  termed  supercritical  flow  (subcritical  flow  occurs  when 
the  flow  is  subsonic  at  all  points  in  the  flowfield;  critical 
flow  occurs  when  M=1  at  one  point  in  the  flowfield) .  Sub¬ 
sonic  flow  accelerating  around  the  airfoil  surface  creates 
regions  where  M>1.  Expansion  waves  leaving  the  airfoil 
surface  will  reflect  off  the  sonic  line  as  compression  waves 
(figure  1).  Some  of  the  compression  waves  will  reach  the 
airfoil  surface  while  others  will  reach  and  strengthen  the 
shock  wave.  The  shock  wave  is  strongest  at  the  airfoil 


surface  and  weakens  to  zero  strength  near  the  top  of  the 
supersonic  flow  region.  Further  complicating  the  flowfield 
is  the  interaction  of  the  shock  wave  and  the  boundary  layer. 
If  the  shock  wave  is  strong  enough  the  boundary  layer  may 
become  separated  from  the  airfoil  surface  thus  extending  the 
viscous  forces  beyond  the  thin  boundary  layer  and  signifi¬ 
cantly  altering  the  inviscid  flowfield.  Separation  of  the 
boundary  layer  at  the  airfoil  trailing  edge  may  also  alter 
the  flowfield.  In  addition  to  creating  a  local  deviation 
from  inviscid  flow  theory,  trailing  edge  separation  can  cause 
a  change  in  circulation  about  the  airfoil  and  thus  a  change 
in  lift.  This  in  turn  can  affect  the  shock  wave  location. 
Finally,  trailing  edge  separation  may  interact  with  a  shock 
wave  and  change  the  pressure  distribution  about  the  airfoil 
surface. 

It  should  be  expected  that  a  phenomenon  as  complicated 
as  the  transonic  airfoil  problem  requires  a  mathematical 
model  equally  complex  if  the  problem  is  to  be  modelled  accu¬ 
rately.  The  governing  equations  for  the  transonic  airfoil 
problem  do  indeed  contain  several  outstanding  mathematical 
complexities  when  the  flowfield  is  supercritical.  First,  if 
the  problem  is  posed  as  a  steady-state  problem  the  governing 
differential  equations  are  of  the  mixed  elliptic-hyperbolic 
type,  being  elliptic  for  subsonic  flow  regions  and  hyperbolic 
for  supersonic  flow  regions.  Second,  the  equation  of  motion 
for  transonic  flow  is  nonlinear.  As  the  freestream  Mach 
number  approaches  unity  the  nonlinear  terms  become  large. 


Thus  near  Mach  one  the  equ«.  i.-.ons  cannot  be  linearized  if  the 
problem  is  to  be  accurately  modelled.  A  third  problem  exists 
when  attempting  to  account  for  the  presence  of  shock  waves. 
Shock  waves  exist  in  flowfields  as  physical  discontinuities 
and  must  be  modelled  as  such.  Including  shock  waves  in  a 
model  requires  the  application  of  the  Rankine-Hugoniot 
relationships  across  the  shock  wave  as  well  as  satisfying  the 
original  governing  equations  throughout  the  rest  of  the 
f lowf ield. 

Two  fundamental  approaches  have  commonly  been  used  to 
analyze  the  transonic  airfoil  problem.  One  is  called  the 
inverse  approach.  Here  the  shock  wave  data  is  specified  and 
the  airfoil  shape  and  details  of  the  flow  within  the  super¬ 
sonic  region  are  unknown.  This  approach  has  the  attractive 
features  of  using  little  conqputer  time  and  requiring 
relatively  little  computer  storage  space.  The  shock  charac¬ 
teristics  are  iterated  in  this  method  until  the  solution  for 
a  specified  airfoil  shape  is  determined.  The  other  approach 
is  called  the  direct  method.  Here  the  airfoil  shape  is  given 
and  all  of  the  details  of  the  flow  within  the  supersonic 
region  are  unknown.  The  direct  method  requires  much  more 
computing  time  and  storage  space  than  the  inverse  method. 

It  stands  to  reason  that  a  set  of  governing  equations 
which  accurately  models  the  transonic  airfoil  problem  will  be 
very  complex.  However,  mathematical  models  which  omit  or 
approximate  some  of  the  phenomena  in  the  f lowf ield  will  be 
less  complex  and  easier  to  work  with.  For  example,  for  flows 


with  extensive  separated  regions  inviscid  flow  theory  must  be 
combined  with  viscous  theory  to  correctly  predict  the  flow- 
field.  This  is  a  complex  model  which  can  be  difficult  to 
analyze.  With  attached  flows  or  flowfields  containing  only 
locally  separated  regions,  inviscid  theory  alone  may  give 
useful  approximations,  thus  allowing  use  of  a  simpler  model. 
Finally,  assuming  shockless  flow  allows  a  simpler  model  to  be 
used  by  requiring  fewer  equations  to  describe  the  flowfield. 
Such  approximations  may  be  made  in  the  early  development  of  a 
solution  technique  to  simplify  the  task  of  proving  the 
validity  of  the  technique.  Once  it  has  been  shown  to  yield 
useful  approximations,  model  refinements  can  be  made  to  have 
the  model  more  accurately  represent  the  flowfield. 

Exact  solutions  to  the  transonic  airfoil  problem  are 
available  for  only  the  simplest  of  cases  (ref  4).  However, 
several  numerical  methods  have  been  developed  which  yield 
approximate  solutions. 

Numerical  Methods  Currently  Available 

Three  types  of  numerical  methods  have  generally  been 
used  to  analyze  the  transonic  airfoil  problem.  One  of  these 
is  the  time-dependent  method.  Here  are  some  of  the  diffi¬ 
culty  in  solving  a  mixed  elliptic-hyperbolic  differential 
equation  is  overcome  by  retaining  the  time  dependent  terms  in 
the  governing  equation.  That  is,  the  problem  is  posed  as  an 
unsteady  one.  This  makes  the  equations  hyperbolic  with 
respect  to  time.  The  original  boundary  value  problem  is 


recast  as  an  initial-value/boundary  value  problem.  A  finite 
difference  scheme  is  then  used  to  integrate  the  equations 
forward  in  time.  If  the  boundary  conditions  are  time- 
independent  or  become  so,  the  solution  will  coverge  to  a 
steady-state  solution  equivalent  to  that  for  the  original 
steady  boundary  value  problem.  Using  the  time-dependent 
approach  assures  that  the  problem  is  properly  posed  math¬ 
ematically.  Magnus  and  Yoshihara  (ref  16)  have  been 
particularly  successful  with  this  method,  having  applied  it 
to  several  airfoil  shapes.  Their  results  were  the  first 
reported  solutions  for  flows  with  imbedded  shock  waves  which 
conqpared  favorably  with  experimental  data.  Grossman  and 
Moretti  (ref  7)  have  also  used  time-dependent  methods  with 
notable  success. 

A  second  type  of  numerical  solution  technique  is  the 
relaxation  methods.  Relaxation  methods  are  a  standard 
procedure  for  solving  boundary  value  problems  described  by 
elliptic  partial  differential  equations.  Relaxation  can  be 
viewed  as  the  recasting  of  a  elliptic  differential  equation 
into  a  parabolic  or  hyperbolic  differential  equation.  Murman 
and  Cole  (ref  20),  Murman  and  Krupp  (ref  21),  Jameson  and 
South  (ref  12),  and  Garabedian  and  Korn  (ref  6)  have,  among 
others,  successfully  applied  relaxation  methods  to  the  tran¬ 
sonic  airfoil  problem.  The  essential  feature  of  their  work 
is  the  accounting  for  the  local  physical  nature  of  the  flow- 
field  while  applying  a  finite  differencing  formula  to  the 
governing  equations  describing  the  flowfield.  Most 
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relaxation  methods  have  been  applied  using  transonic  small 
disturbance  theory.  Steger  and  Lomax  (ref  22)  have  applied 
relaxation  using  both  small  disturbance  theory  and  full 
inviscid  theory. 

Finally,  various  approximation  methods  have  been  applied 
to  the  transonic  airfoil  problem.  These  methods  employ  math¬ 
ematical  simplifications  to  reduce  partial  differential 
equations  to  ordinary  differential  equations.  The  approach 
developed  in  this  study  is  of  this  type.  Approximations 
methods  have  been  successfully  applied  using  both  the  inte¬ 
gral  method  and  the  method  of  integral  relations. 

The  integral  method  recasts  the  governing  equations  into 
a  set  of  integral  relations  through  the  u3e  of  Green's 
theorem.  An  approximation  is  then  used  to  represent  the 
variation  of  a  dependent  variable  in  one  direction.  The 
equation  is  then  integrated  in  that  direction  leaving  an 
ordinary  differential  equation  with  respect  to  the  other 
independent  variable.  This  method  has  not  been  found  to 
consistently  produce  results  comparing  favorably  with 
experimental  data. 

The  method  of  integral  relations  differs  from  the  above 
method  in  that  the  flowfield  is  divided  into  strips,  thus 
discretizing  one  of  the  independent  variables.  Separate 
approximations  for  the  dependence  of  the  dependent  variables 
with  respect  to  the  discretized  variable  are  then  used  for 
each  strip.  A  matching  condition  is  imposed  at  the  strip 
boundaries.  This  method  has  been  successfully  used  by  Tai 
(ref  23),  Melnik  and  Ives  (ref  18),  Holt  and  Masson  (ref  9), 
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Crown  (ref  2)  and  Dorodnitsyn  (ref  3).  Informative  summaries 
of  numerical  work  on  the  transonic  airfoil  problem  have  been 
written  by  Murman  (ref  19),  Taylor  (ref  24),  and  Yoshihara 
(ref  25). 

Statement  of  Problem 

The  objective  of  this  study  is  to  develop  and  evaluate  a 
simple  numerical  solution  technique  applicable  in  both  the 
subcritical  and  supercritical  regimes.  The  method  should 
provide  a  good  approximation  to  the  flow  field  being  modelled 
and  should  use  relatively  little  computer  time.  The  method 
will  first  be  applied  to  the  simple  case  of  incompressible 
flow  around  the  unit  circle  airfoil.  The  method  will  then  be 
tested  on  compressible  flow  around  the  unit  circle  airfoil. 
Through  the  application  of  an  appropriate  transformation  the 
circular  airfoil  results  will  be  applied  to  Joukowski  airfoil 
shapes  for  the  incompressible  flowfield  case.  To  reduce  the 
complexity  of  the  problem  only  shockless  flowfields  will  be 
studied.  Flowfields  will  be  assumed  to  be  inviscid  and  all 
flows  will  be  irrotational  and  isentropic. 

Approach  to  the  Problem 

The  full  inviscid  potential  theory  and  the  irrotational- 
ity  condition  will  provide  the  set  of  governing  equations  to 
be  used  in  this  study.  The  boundary  conditions  to  be  applied 
to  the  problem  are  no  flow  normal  to  the  airfoil  at  the  air¬ 
foil  surface  and  uniform  flow  at  infinity.  A  major  facet  of 
the  problem  formulation  will  be  a  coordinate  transformation 


applied  to  the  governing  equations  and  boundary  conditions. 
The  new  coordinate  system  will  simplify  the  geometry  of  the 
problem  and  will  make  application  of  the  infinity  boundary 
condition  easier.  Of  the  three  types  of  numerical  solution 
techniques  mentioned  it  was  decided  that  one  of  the  approxi¬ 
mation  methods  would  be  used  in  this  study.  This  choice  was 
made  based  upon  the  desire  to  keep  the  numerical  model 
relatively  simple.  It  has  been  shown  by  McCracken  (ref  17) 
that  useful  results  can  be  obtained  when  approximate  methods 
are  used  in  a  simplified  system  model.  Specifically  the 
method  of  lines,  a  version  of  the  method  of  integral  rela¬ 
tions,  will  be  used  in  this  study.  This  solution  technique 
has  been  used  by,  among  others,  McCracken  (ref  17),  Klunker, 
Davis  and  South  (ref  15),  Hamilton  (ref  8),  and  Melnik  and 
Ives  (ref  18).  The  method  of  lines  will  be  used  to  reduce 
the  set  of  governing  partial  differential  equations  to  a  set 
of  ordinary  first  order  differential  equations.  This  new  set 
will  then  be  numerically  integrated  to  yield  a  velocity  dis¬ 
tribution  along  the  airfoil  surface.  The  solution  technique 
will  be  evaluated  by  comparing  the  computed  velocity  distri¬ 
butions  with  published  data. 

The  unit  circle  airfoil  was  chosen  because  the  exact 
solution  for  incompressible  flow  around  it  is  known.  It  thus 
presented  itself  as  one  way  in  which  the  accuracy  of  the 
solution  technique  developed  in  the  study  could  be  tested. 
Another  reason  the  unit  circle  airfoil  was  chosen  is  that  the 
Joukowski  transformation  may  be  applied  to  it  thus  allowing 
airfoil-like  shapes  to  be  studied. 


II .  Development  of  the  Solution  Technique 

There  are  several  methods  available  with  which  to 
analyze  the  2-D  transonic  airfoil  problem.  Regardless  of  the 
method  used  certain  choices  must  be  made  which  will  define 
the  approach  used  to  work  the  problem.  First  a  set  of 
governing  equations  mast  be  chosen.  Boundary  conditions  must 
also  be  specified.  Second,  a  coordinate  system  must  be 
chosen  within  which  to  work  the  problem.  Third,  a  numerical 
solution  technique  must  be  developed  to  solve  the  set  of 
governing  equations. 

When  choosing  a  set  of  governing  equations  a  trade-off 
is  made  between  the  accuracy  desired  in  the  system  model  and 
the  level  of  coiqplexity  acceptable.  For  inviscid  flow 
analysis  one  of  three  sets  of  relationships  have  commonly 
been  used  to  mathematically  model  the  2-D  transonic  airfoil 
problem.  They  differ  primarily  in  hew  each  accounts  for  the 
presence  of  shock  waves. 

One  set  consists  of  the  isentropic  relationship 


along  with  the  unsteady  form  of  the  conservation  equations 
for  mass  and  momentum 


+  (VpQ)  «  0  (2) 

If  +  (Q-V)  Q  -  -7  (fj 

These  are  applied  throughout  the  flowfield  except  across  the 
shock,  where  the  Rankine-Hugoniot  relationships  are  used. 


The  number  of  equations  in  this  formulation  however  makes  it 
difficult  to  work  with.  Such  a  set  of  equations  has  been 
used  by  Grossman  and  Moretti  (ref  7). 

A  second  set  of  equations  consists  of  the  steady-state 
form  of  the  previous  set  written  in  divergence  form.  This 
set  can  be  applied  across  a  shock  wave  thus  avoiding  the  need 
for  special  relationships  across  the  shock  wave.  However, 
this  set  of  equations  is  of  the  mixed  type.  This  set  has 
been  employed  in  2-D  transonic  airfoil  analysis  by  Tai  (ref 
23) . 

A  third  set  of  equations  may  be  fomulated  starting  with 
the  assumption  that  a  good  approximation  to  the  flowfield  can 
be  obtained  without  including  shock  waves  in  the  model.  This 
assumption  is  compatible  with  the  requirement  for  the  flow  to 
be  isentropic  since  isentropic  flow  implies  either  weak 
shocks  or  no  shocks  present  in  the  flowfield.  Based  upon  the 
desire  to  start  with  as  simple  a  model  as  possible  it  was 
decided  to  use  this  set  and  neglect  shock  waves.  Therefore 
the  equations  used  in  this  study  model  the  2-D  transonic 
airfoil  problem  assuming  steady,  inviscid,  irrotational, 
isentropic  flow.  Given  these  assumptions,  inviscid  flow 
theory  and  the  irrotationality  condition  were  used  to  define 
a  set  of  governing  equations  in  cylindrical  coordinates.  The 
nondimens ionalized  equations  are  stated  here.  The  interested 
reader  is  referred  to  the  Appendix  for  the  details  of  their 
development.  The  irrotationality  equation  is 


and  the  inviscid  full  potential  equation  of  motion  is 


3 v  2  Euv  rav  _  B-Du*-Av*  r au 

3 e  *  B-Au*-Dv*  3 r  ~  B-Au*-Dv*  3r 


For  flow  around  the  unit  circle  airfoil  the  boundary  condi¬ 
tions  were  that  the  normal  velocity  component  at  the  airfoil 
surface  equal  zero  and  that  uniform  flow  conditions  exist  at 
infinity.  In  nondimensionalized  form  the  zero  velocity 
boundary  condition  was  written 

u(l , 9)  *  0  (6) 

and  the  infinity  boundary  condition  was  written 

u*(-,e)  +  v*(«,e)  »  l  (7) 

The  decision  to  use  the  unit  circle  airfoil  made  the 
choice  of  a  cylindrical  coordinate  system  obvious.  It  was 
noted  though,  that  this  choice  could  be  improved  by  applying 
the  simple  inversion  coordinate  transformation 


n 


(8) 


The  effect  of  this  transformation  is  shown  in  figure  2.  The 
advantages  gained  by  its  application  are  two-fold.  First, 
the  domain  of  the  problem  is  reduced  to  a  finite  region. 
Second,  the  transformation  solves  the  problem  of  where  to 
apply  the  infinity  boundary  condition.  The  application  of 
this  coordinate  transformation  to  the  governing  equations  and 
boundary  conditions  is  shown  in  the  Appendix.  The  results 
are  that  the  irrotationality  condition  becomes 


and  the  equation  of  motion  becomes 


3V  —  Av1tDu1-B  n3U  +  2Euv  na v  (10) 

30  '  Au*+Dv*-B  3,1  Au*+Dv*-B  3n 

The  boundary  condition  of  no  normal  velocity  component  at  the 
airfoil  surface  remains  unchanged  by  the  transformation.  The 
infinity  boundary  condition  becomes 

u* (0, e)  +  v*(0,e)  *  1  (11) 

There  remained  the  choice  of  which  numerical  solution 
technique  to  use  to  solve  this  set  of  equations.  It  was 
decided  to  use  the  method  of  lines  because  of  its  ease  of 
application  to  the  governing  equations  used  in  this  study. 

The  method  of  lines  is  used  to  reduce  the  partial 
differential  equations  to  first  order  ordinary  differential 
equations.  This  is  done  by  allowing  the  derivatives  with 
respect  to  n  to  be  represented  by  a  finite  difference 
approximation.  As  illustrated  in  Figure  3  the  n  domain  is 
divided  by  N  concentric  lines  equally  spaced  between  n  =*  0 
and  n  *  1  such  that 

”i  =  I?!'  1  *  2 ,  3 ,  ...»  N,  N+l  (12a) 

with  ni  »  0  and  nN+2  “  1  (12b) 

Along  the  N  lines  the  derivatives  with  respect  to  n  are 
approximated  by  a  central  difference  formula.  Details  of  the 
formula  and  how  it  was  used  in  the  study  are  shown  in  the 
study  are  shown  in  the  Appendix.  The  result  is  that  along 
the  itl1  line  the  governing  equations  are  written 


.2 
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v  -  1=1 

vi  N+l 


vi+l~vi-l 


2 
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(14) 


Writing  these  equations  for  each  of  the  N  lines  results 
in  a  set  of  2N  simultaneous  equations  involving  2N+4  unknowns, 
u^  and  v^,  i=*l,  2,  . . . ,  N,  N+l,  N+2 .  Clearly  four  more 
relationships  are  needed  to  make  the  problem  determinate. 

Three  of  these  are  obtained  from  the  boundary  conditions  as 


u  (0 , 0 )  = 

-cos  0 

(15) 

v(O,0)  * 

sin  0 

(16) 

u(l ,0 )  = 

0 

(17) 

The  fourth  relationship  is  derived  by  applying  the  irrotation- 
ality  condition  along  the  airfoil  surface.  Here  it  is  known 
that  the  radial  velocity  component  u  is  equal  to  zero  for  all 

9  U 

e.  This  implies  that  along  the  airfoil  surface  —  is  equal  to 

a  0 

zero  for  all  0.  Thus,  for  n“l  (i.e.,  along  the  airfoil 
surface)  one  can  write 


a  u 
90 


0 


(18) 


After  substituting  for  n  and  applying  a  backward  difference 
formula  to  approximate  the  derivative  with  respect  to  n , 
equation  (18)  is  rewritten 
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To  summarize,  the  problem  now  consists  of  2N  first  order 
nonlinear  ordinary  differential  equations  and  one  algebraic 
relationship  along  with  three  boundary  conditions,  two  of 
these  applied  at  n  *  0  for  all  s,  and  the  third  applied  at 
n  *  1  for  all  9.  All  that  remained  was  to  develop  a  scheme 
with  which  to  integrate  these  equations  and  to  state  the 
initial  conditions  to  be  used. 

For  the  sake  of  simplicity  it  was  decided  to  use  a 
"canned"  integration  package  called  DVERK  which  was  available 
in  the  International  Mathematical  and  Scientific  Library 
( IMSL )  of  computer  programs.  DVERK  finds  approximations  to 
the  solution  of  a  system  of  first  order  ordinary  differential 
equations  with  initial  conditions.  The  program  uses  Runge- 
Kutta  formulas  of  order  5  and  6.  As  applied  in  this  study, 
DVERK,  given  a  set  of  initial  conditions,  would  integrate  the 
2N  ordinary  differential  equations  with  respect  to  9  by  some 
amount  A9.  The  integration  stepsize  could  be  chosen  by  the 
user  to  be  whatever  size  was  necessary  to  obtain  satisfactory 
results.  In  this  study  A9  was  typically  one  degree.  The 
product  of  each  integration  step  consisted  of  the  normal  and 
tangential  velocity  components  on  each  line  at  the  angle  the 
problem  has  been  integrated  to.  The  algebraic  relationship 
developed  by  applying  the  irrotationality  condition  at  the 
airfoil  surface  was  then  used  to  compute  the  velocity  at  the 
airfoil  surface  at  that  same  angle.  The  velocities  on  each 


line  as  well  as  at  the  airfoil  surface  then  served  as  the 
initial  conditions  for  DVERK  to  integrate  another  amount  A9. 
By  integrating  with  repsect  to  9  from  zero  to  ninety  degrees 
a  velocity  profile  along  the  airfoil  surface  would  be 
generated. 

Next,  a  scheme  for  providing  a  set  of  correct  initial 
conditions  at  9  *  0°  had  to  be  developed.  This  was  not 
necessary  for  the  incompressible  flowfield  case  since  exact 
solutions  for  this  case  are  available.  However,  a  scheme  had 
to  be  developed  which  would  provide  a  correct  set  of  initial 
conditions  at  9  -  0°  for  compressible  flowfield  cases.  The 
following  discussion  outlines  the  scheme  developed  to  obtain 
these  initial  conditions  for  compressible  flowfield  cases. 

Correct  initial  conditions  were  obtained  by  starting 
with  an  estimated  set  of  initial  conditions  at  9  »  O',  a 
frees tream  Mach  number  slightly  greater  than  zero  (typically 
0.1),  and  using  from  three  to  seven  lines.  Using  this  set  of 
estimated  initial  conditions  the  governing  equations  were 
integrated  to  9  *  90°.  After  marching  ninety  degrees,  normal 
and  tangential  velocity  components  were  available  on  each 
line  at  9  =  90°  as  well  as  at  the  airfoil  surface.  For  the 
flowfields  analyzed  in  this  study  the  normal  velocity  com¬ 
ponents  are  zero  at  9  =*  90®.  This  ideal  velocity  profile  was 
never  realized  however  after  having  numerically  integrated 
one  time  to  9  =  90®  when  starting  with  an  estimated  set  of 
initial  conditions  at  9  *  0®.  It  was  assumed  however  that 
there  did  exist  some  velocity  distribution  at  9  *  0®  which 
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if  used  as  a  set  of  initial  conditions,  would  allow  the 
integration  to  yield  a  correct  velocity  distribution  (i.e. 
one  matching  the  ideal)  at  9  =  90®.  The  problem  was  to  find 
this  velocity  profile. 

It  was  decided  that  this  could  best  be  accomplished  by 
systematically  adjusting  the  normal  velocities  at  e  =0°  to 
obtain  normal  velocities  of  zero  at  0  =  90®.  A  first  order 
Newton-Raphson  iterative  scheme  was  used  to  determine  the 
changes  required.  Using  this  scheme,  each  normal  velocity 
component  at  0  =  0®  was,  in  turn,  varied  a  small  amount  from 
its  original  estimated  value.  DVERK  was  then  used  to  inte¬ 
grate  the  governing  equations  through  ninety  degrees.  The 
change  in  each  normal  velocity  component  at  0  =  90®  as  a 
result  of  varying  each  normal  velocity  component  at  0  *  0® 
was  thus  calculated.  These  changes  were  used  to  compute 
rates  of  change  of  velocities  at  8  =  90®  as  functions  of 
small  velocity  changes  at  0  *  0®.  The  rates  were  then  used 
in  an  iterative  Newton-Raphson  scheme  to  determine  the  amount 
each  normal  velocity  component  at  0  =  0®  had  to  be  varied 
from  its  original  estimated  value  to  produce  normal  velocity 
components  of  value  zero  at  0  =  90®. 

If  the  problem  of  computing  a  correct  normal  velocity 
profile  at  0  a  0®  were  linear,  only  one  iteration  of  the 
Newton-Raphson  scheme  would  be  needed.  However,  the  problem 
is  non-linear  and  several  iterations  of  the  Newton-Raphson 
scheme  were  needed.  Use  of  the  Newton-Raphson  scheme  imposed 
one  constraint  upon  the  original  estimate  of  the  normal 


velocity  profile  at  0  =0°.  Since  it  is  a  linear  method 
being  used  to  solve  a  non-linear  problem,  the  original 
estimate  could  not  be  too  far  from  the  final  values.  This 
constraint  was  overcome  by  simply  insuring  that  only  small 
changes  were  needed  to  the  original  estimate.  This  was 
accomplished  by  starting  with  the  known  exact  velocity  pro¬ 
file  for  incompressible  flow  and  solving  a  flowfield  problem 
at  only  a  slightly  higher  Mach  number.  Then  using  this  new 
profile  the  Mach  number  was  increased  slightly  again  and  the 
problem  solved  again.  This  incremental  process  was  repeated 
until  a  solution  was  obtained  for  the  desired  freestream  Mach 
number . 

In  this  study,  a  flowfield  was  analyzed  starting  with 
the  incompressible  velocity  profile  at  0  *  0°  and  solving  for 
the  velocity  profile  (and  correct  initial  conditions)  at.  M  ■ 
0.1.  Results  for  M  =  0.1  were  then  used  as  a  starting  point 
in  analyzing  the  flowfield  at  M  =  0.2,  and  so  forth  up  to  the 
desired  Mach  number.  This  scheme  proved  to  be  convergent  for 
subcritical,  critical  and  slightly  supercritical  flowfields. 
For  these  cases  the  normal  velocities  at  0  =  90°  could  be 
made  to  approach  zero  using  as  few  as  four  lines. 


III.  Results 


The  results  and  conclusions  of  the  study  are  presented 
in  order  of  increasing  complexity  of  the  flowfield  being 
analyzed.  The  results  for  the  case  of  the  unit  circle 
airfoil  in  an  incompressible  flowfield  are  first  presented 
followed  by  the  results  for  a  symmetrical  Joukowski  airfoil 
in  an  incompressible  flowfield.  Results  for  the  unit  circle 
airfoil  in  a  compressible  flowfield  are  then  presented. 

Incompressible  Flowfield  Results 

Incompressible  flow  about  the  unit  circle  airfoil  was 
chosen  as  the  first  application  of  the  solution  technique  for 
several  reasons.  First,  this  simple  case  has  a  known  exact 
solution  (ref  14).  The  results  from  the  solution  technique 
developed  in  this  study  could  be  compared  to  this  exact 
solution.  This  comparison  would  be  useful  as  an  indication 
of  the  accuracy  of  the  solution  technique  developed.  A 
second  reason  lay  in  the  fact  that  the  governing  equations 
for  this  case  were  known  to  have  a  high  degree  of  numerical 
stability  due  to  their  linearity  in  the  variables  u  and  v. 

The  irrotationality  condition  is  always  linear  in  these 
variables.  However,  the  equation  of  motion  is  linear  only 
for  the  case  of  incompressible  flow.  For  this  case  the 
equations  are 


Starting  with  the  boundary  conditions  as  previously 
stated  the  solution  technique  was  numerically  intergrated 
around  the  airfoil  from  0  =  0  to  0  =*  90*.  Because  of  the 
symmetry  of  the  problem  all  of  the  flowfield  was  known  if  the 
flowfield  in  this  quandrant  was  known.  The  numerical  inte¬ 
gration  generated  a  velocity  distribution  along  the  airfoil 
surface.  This  data  was  then  compared  to  the  exact  solution. 
Figure  4  shows  the  velocity  distribution  at  the  surface  of 
the  unit  circle  airfoil  in  an  incompressible  flowfield  along 
with  the  exact  solution.  Three  lines  were  used  for  this 
case.  It  was  found  that  for  incompressible  flow  as  few  as 
three  lines  were  needed  for  accurate  results.  Increasing  the 
number  of  lines  did  not  increase  the  accuracy  when  analyzing 
incompressible  flowfields. 

Having  solved  the  case  of  incompressible  flow  about  the 
unit  circle  airfoil,  the  Joukowski  transformation  was  applied 
to  that  flowfield  to  analyze  the  case  of  incompressible  flow 
about  a  Joukowski  airfoil.  Figures  5  and  6  show  the  velocity 
distribution  for  two  typical  Joukowski  airfoil  profiles  as 
obtained  by  the  solution  technique  developed  in  this  study. 

Compressible  Flowfield  Results 

Unit  circle  airfoil  surface  velocity  profiles  were 


obtained  over  a  range  of  Mach  numbers.  Figure  7  shows  the 
results  for  a  subcritical  flowfield  case.  Figure  8  shows  the 


results  for  the  critical  flowfield  case  where  M  =  0.42. 


Compressible  flowfields  with  freestream  Mach  numbers  up 
to  0.45  were  successfully  studied.  Attempts  to  analyze  flow- 
fields  with  freestream  Mach  numbers  higher  than  0.45  failed 
however.  Neither  the  Newton-Raphson  linear  scheme  nor  a  non¬ 
linear  least  squares  scheme  (ref  5)  would  produce  convergent 
results.  Varying  the  numerical  integration  step  size  and  the 
number  of  lines  did  not  solve  the  problem.  According  to 
previously  published  reports  in  which  the  method  of  lines  was 
used  (ref  13),  the  difficulty  is  due  to  the  inherent  insta¬ 
bility  present  in  the  governing  equations  when  the  method  of 
lines  is  used  with  a  central  difference  approximation  for  the 
derivatives.  Central  differencing  is  appropriate  for  approx¬ 
imating  derivatives  in  elliptic  equations  (subcritical  flow) 
where  this  differencing  can  yield  convergent  results.  But 
for  hyperbolic  differential  equations  (supercritical  flow) 
central  difference  approximations  are  unstable  and  will  not 
yield  convergent  results. 

Figure  9  shows  the  velocity  distribution  for  the  highest 
freestream  Mach  number  that  the  solution  technique  would 
successfully  analyze.  At  this  Mach  number,  M.  *  0.45,  it 
should  be  noted  that  the  maximum  local  Mach  number  is  in  the 
flowfield  is  1.45.  Therefore  the  solution  technique  can  be 
said  to  be  useful  in  analyzing  supercritical  flowfields  to  at 
least  a  limited  extent. 
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concluded  that  the  solution  technique  produces  accurate 
results  when  analyzing  subcritical  flow  around  circular 
airfoil  shapes.  The  method  also  gives  good  results  for 
slightly  supercritical  flowfields.  However,  when  analyzing 
supercritical  flowfields  around  circular  shaped  airfoils  the 
method  breaks  down  for  M«>.45.  It  is  concluded  that  the 
finite  differencing  scheme  used  in  this  study,  i.e.  the 
method  of  lines  with  central  differencing,  is  appropriate 
when  the  system  governing  equations  are  of  the  elliptic  type 
This  occurs  for  subcritical  flows  around  circular  airfoil 
shapes.  The  scheme  is  however  incorrect  when  analyzing 
supercritical  flowfields  because  central  differencing  alone 
cannot  account  for  the  local  nature  of  the  flowfield  in 
regions  of  supercritical  flow. 


Recommendations 

The  following  recommendations  are  proposed  for  any 
future  work  in  modelling  flowfields  with  the  method  of 
lines: 

1.  A  finite  differencing  and  integration  scheme  using 
the  method  of  lines  should  be  developed  which  would  account 
for  the  local  physical  nature  of  flowfields  in  the  transonic 


regime.  The  integration  scheme  should  be  capable  of  sensing 
when  the  flowfield  is  transitioning  from  subsonic  to  super¬ 
sonic  flow.  The  scheme  should  then  be  capable  of  applying 
the  appropriate  differencing  formulas  in  response  to  the  type 
of  flowfield  sensed. 

2.  An  inversion  transformation  such  as  the  type  used 
in  this  study  should  be  developed  which  will  allow  the  solu¬ 
tion  technique  used  to  be  applied  to  a  variety  of  airfoil 
shapes.  The  benefits  of  such  a  transformation  include  easier 
application  of  the  flowfield  boundary  conditions  as  well  as 
extending  the  utility  of  the  technique. 
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Figure  4 


Solution  Obtained  for  the  Velocity  Distribution  at  the 
Surface  of  the  Unit  Circle  for  M^O  as  Compared  to  the 
Exact  Solution.  N=3 
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Figure  5 


Solution  Obtained  for  the  Velocity  Distribution  at  the 
Surface  of  a  2%  Thick  Symmetrical  Joukowski  Airfoil  for 
M.=0  as  Compared  to  the  Exact  Solution.  N=3 
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Figure  6 


Solution  Obtained  for  the  Velocity  Distribution  at  the 
Surface  of  a  14%  Thick  Symmetrical  Joukowski  Airfoil  for 
M«»0  as  Conqpared  to  the  Exact  Solution.  N«3 


\ 


31 


Solution  Obtained  for  the  Velocity  Distribution  at  the 
Surface  of  the  Unit  Circle  for  the  Sub-Critical  Plowfield 
Case  of  M .”0.4  as  Compared  to  Published  Date.  N®»4 
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Figure  3 


Solution  Obtained  for  the  Velocity  Distribution  at  the 
Surface  of  the  Unit  Circle  for  the  Critical  Flowfield  Case 
of  M^0.42  as  Compared  to  Published  Data.  N*4 


This  appendix  provides  a  detailed  account  of  the 
development  of  the  solution  technique  used  in  this  study. 

The  governing  equations  and  boundary  conditions  are  first 
developed.  Then  the  equations  are  nondimensionalized  and 
arranged  into  a  more  concise  form.  Next  a  coordinate 
transformation  is  applied  to  the  nondimensionalized  governing 
equations  and  boundary  conditions.  Finally,  the  method  lines 
is  used  to  reduce  the  governing  equations  from  partial  to 
ordinary  differential  equations. 


Development  of  the  Governing  Equations 

The  set  of  governing  equations  used  in  this  study  are 
used  to  model  steady,  inviscid,  irrotational,  isentropic 
flow.  As  a  starting  point  the  continuity  equation 

V  •  ( pQ)  -  0  (22) 


the  equation  of  motion 

(Q  •  v )  Q  »  -v  (  E  j 

and  the  equation  of  state 


(23) 


p  ■  kpT  (24) 

are  used  from  which  the  final  equations  will  be  derived.  The 


speed  of  sound,  c,  given  by  c* 


LE 

»p 


,  may  be  introduced 


35 


into  equation  (23)  which,  when  then  combined  with  equation 


(22),  gives 


Q  •  7 


(§’)- 


c1  7  •  Q  *  0 


letting  Q* 


Q  •  Q. 


The  definition  of  c*  may  again  be  used  to  combine  equations 
(23)  and  (24). 

The  resulting  equation  can  be  integrated  once  to  give 

^r  +  f  -  ik  +  T-  <26> 

where  the  constant  of  integration  is  chosen  to  match  the 
conditions  of  undisturbed  flow,  i.e.,  Q  =  (-Q-  cos  e, 

Q„  sin  e)  and  c  *  c„. 

With  the  assumption  of  irrotationality,  a  potential 
function  »  may  be  introduced  into  equations  (25)  and  (26) 


with 


—  »  Q  a 

ar  r 

i  14.  m  Q 
r  ae  e 


-Q«  cos  9 


Q»  sin  e 


(27a) 


(27b) 


Equation  (25)  may  then  be  rewritten  as 

l  /  a*  3Q*  .  l  a*  aQ* 


ar  ar  r*  ae  ae 


with 


fi±y  .  1  (9±\* 

\ar /  *  r*  e/ 


Equation  (26)  may  also  be  rewritten  as 
c*  -  c*  -  (Q1  -  Qi) 


Equation  (28)  is  the  equation  of  motion  used  for  compressible 
potential  flow.  The  assunqption  of  irrotational  flow  allows 
the  use  of  a  potential  function.  The  irrotationality 
condition  is  expressed  as 
7  x  Q  ■  0 
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(31) 
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One  boundary  condition  to  be  used  is  the  assumption  of 
uniform  flow  at  infinity.  This  boundary  condition  will  be 
applied  in  the  form 

Q*  a  (-Q„  cos  9)*  +  (Q^  sin  9)* 

=  Q*  +  Q*  =  U*  +  V*  (32) 

r  0 

The  second  boundary  condition  is  that  there  be  no  normal 
velocity  component  at  the  airfoil  surface,  i.e. 

Ui  «  0  (33) 

along  the  airfoil  surface. 


Nondlmensionalization 

Defining  ♦  “  o 

the  equation  of  motion 
as 


(i  _m«  (q*-l)j 

with  q*  » 

ar  r  * 


,  q  * 
may  be 


§  and  M.  -  (34) 

written  in  nondime ns ional  form 


mi  /  . 

-  n  iai  +  I  ii  iaM 

2  \ar  ar  r*  ae  99  / 

(35) 

(36) 


and 


c* 

Qi 


1 

Mi 


(■ 


M- 


<q-  -1)1 


(37) 


Equation  (37)  when  expanded  can  be  written 


i-ii.  +  I  H±+k±±)  , 

ar*  r*39*  r  3r  I 


Mi 


J-il  +  1  11  ii  ?*♦  .1  11  (14.)* 
ar*  r*  ae  arae  r»  3r  ^39/ 


1 

r* 


111 

ae* 


+ 

(38) 
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Assuming  potential  flow  allows  one  more  simplification  to  be 
made  by  allowing 

If  -  u  and  iff  -  v  (39) 

with  u  being  the  nondimens ional  radial  velocity  component  and 
v  being  the  nondimens ional  tangential  velocity  component. 
Applying  these  definitions  to  equation  (38)  yields,  after 
some  rearrangement. 


Finally,  this  equation  may  be  rearranged  into  the  form 


9  v 

Av*+Du*-B  au 

2Euv  3u 

(41) 

30 

B-Dv* -Au*  r  3r 

B-Dv*  -Au*  3  0  U 

with 

A  = 

Mi 

(y-1) 

2 

(42a) 

B  = 

M* 

<^>  +  i 

(42b) 

D  * 

Mi 

(xr'  *  »i 

and 

(42c) 

E  * 

Mi 

(42d) 

Equation  (41)  is  the  equation  of  motion  in  the  form  used  in 
this  study.  The  irrotationality  condition  may  be  written  in 


similar  form  as 


I- 


£ 


h 


Nondimensionalizing  the  infinity  boundary  condition, 

Qi  (-Q„  cose)1  (Q.  sin  8)* 


equation  (32),  yields 

U*  V* 

Ql  +  Qi 


Qi 

u*  +  v* 


Qi 


Qi 


(44) 


The  nondimens ional  surface  boundary  condition  simply  becomes 


0 

—  =  u  =*  0 


(45) 


Application  Of  The  Coordinate  Transformation 

To  simplify  the  geometry  of  the  problem  and  the 
application  of  the  boundary  conditions  the  coordinate 
transformation  n  =*  -p  is  employed.  First,  the  transformation 

laws  are  developed.  Assume  that  some  function  f(r,e)  ** 
F(n,8)  exists,  who?^  value  doesn't  change  with  a  change  in 
coordinate  systems.  The  following  equations  can  then  be 
applied: 


3f 

3F  3  f| 

4, 

3F 

38 

(46) 

3r 

3n  3r 

* 

38 

3r 

3  f 

8F  3  n 

+ 

3F 

3  8 

(47) 

38 

3n  38 

38 

38 

3n 

3r 

-i 
*  r7 

(48a) 

3n 

38 

38 

m  — 

3r 

- 

0, 

and 

(48b) 

38 

38 

»  1 

(48c) 

Substituting  equations  (48)  into  equations  (46)  and  (47) 

-  «*  —  (49) 


.  ,  .  3f  -1  3F 

yields  —  -  p-  — 


3F 

3n 
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Equations  (49)  and  (50)  are  the  transformation  laws  needed. 
Applying  them  to  equations  (41)  and  (43)  gives 

a v  _  Av*+Du*-B  au  2Euy  n  av  .  _ . » 

39  =  Au*+Dv*-B  3n  Au*+Dv*-B  an  “  u  '  ; 


Application  of  the  Method  of  Lines 

By  using  the  method  of  lines  the  derivatives  with 
respect  to  n  are  represented  by  a  finite  difference 
approximation.  Hornbeck  (ref  10)  shows  how  the  finite 
difference  approximations  are  derived  as  well  as  the  error 
terms  implicit  in  their  use. 

After  application  of  the  coordinate  transformation  the 
new  domain  is  divided  into  N+l  annular  strips  by  N  lines 
equally  spaced  An  =  apart.  The  value  of  n  along  any  line 


is  then  n^ 


2,  3,  ...,  N,  N+l 


i*l  corresponds  to  n*0 ,  i»N+2  corresponds  to  n«l •  Along  any 
of  the  N  lines  the  first  order  partial  derivative  with 
respect  to  n  of  a  function  F(n»8)  is  approximated  by  the 
second  order  central  difference  formula 


F(n+An>8)  -  F(n-An>9) 
2  An 


i 
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Along  the  unit  circle  n*l •  Here,  since  only  values  of  n  less 
than  1  are  available,  the  partial  with  respect  to  n  is 
approximated  by  the  second  order  backward  difference  formula 


3F(1,3)  m  3F (1 , 6 )  -  4F(1-An,3)  +  F(l-2An,9) 
3  n  2  An 


(55) 


Using  subscripts  corresponding  to  the  value  of  desired, 
the  finite  difference  formulae  may  be  written 


3F 

3n 


,A  -  ™  (Fi+i  -  pi-i) 


for  the  central  difference  approximation  and 


3F 

an 


N+l  /. 

al  "  2  V 


3FN+2  “  4FN+1 


♦  Fh) 


(56) 


(57) 


for  the  backward  difference  approximation. 

Finally,  applying  equations  (56)  and  (57)  to  the  governing 
equations  (51)  and  (52)  yields 


3  v. 

/Av?+Du?-B\ 

/  V 

i 

I  1  1  1 

i-1 

fui+i  "  ui-i) 
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\Auf+Dv£-B  J 

2 

l 11  1  JJ 

(au’«w>-b)  (Vl+1  “  Vl'lj 

for  the  equation  of  motion  and 

.  v  .  <izi>  (v  _  v  ) 
33  Vi  2  \  i+1  Vi-1 J 


-u. 


(58) 


(59) 


for  the  irrotationality  condition 
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